Overcoming the critical slowing down of flat-histogram Monte Carlo simulations: 
Cluster updates and optimized broad-histogram ensembles 
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We study the performance of Monte Carlo simulations that sample a broad histogram in en- 
ergy by determining the mean first-passage time to span the entire energy space of d-dimensional 
ferromagnetic Ising/Potts models. We first show that flat-histogram Monte Carlo methods with 
single-spin flip updates such as the Wang-Landau algorithm or the multicanonical method perform 
sub-optimally in comparison to an unbiased Markovian random walk in energy space. For the 
d — 1, 2,3 Ising model, the mean flrst-passage time r scales with the number of spins N — L'^ as 
T (X N'^L^. The critical exponent z is found to decrease as the dimensionality d is increased. In 
the mean-field limit of infinite dimensions we find that z vanishes up to logarithmic corrections. 
We then demonstrate how the slowdown characterized by z > for finite d can be overcome by 
two complementary approaches - cluster dynamics in connection with Wang-Landau sampling and 
the recently developed ensemble optimization technique. Both approaches are found to improve the 
random walk in energy space so that r oc A*'^ up to logarithmic corrections for the d = 1,2 Ising 
model. 

PACS numbers: 02.70.Rr, 75.10.Hk, 64.60.Cn 



o 



> 
o 

O 



I 

O 

o 



X 



I. INTRODUCTION 

Recently, Wang and Landau P, 01 introduced a Monte 
Carlo (MC) algorithm that simulates a biased random 
walk in energy space, systematically estimates the den- 
sity of states, g{E), and iteratively samples a flat his- 
togram in energy. The bias depends on the total energy 
i5 of a configuration and is defined by a statistical ensem- 
ble with weights w{E) — l/g{E). The idea is that the 
probability distribution of the energy, p{E) — w{E)g{E), 
will eventually become constant, producing a flat energy 
histogram. The algorithm has been successfully applied 
to a wide range of systems 0. 

One measure of the performance of the Wang-Landau 
algorithm and other broad-histogram algorithms is the 
mean first-passage time t, which we define to be the 
number of MC steps per spin for the system to go from 
the configuration of lowest energy to the configuration of 
highest energy. This time has been called the tunneling 
time or round-trip time in earlier work ^ Q . This time 
is the relevant scale for sampling the entire energy space. 
Because the number of possible energy values in the Ising 
model scales linearly with the number of spins N = L'^, 
where d is the spatial dimension, we might expect from 
a simple random walk argument that r oc iV^. However, 
it was recently shown that for flat-histogram algorithms 



T scales as 
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where the exponent z is a measure of the deviation of the 
random walk from unbiased Markovian behavior 4]. In 
analogy to the dynamical critical exponent for canonical 
simulations, the exponent z quantifies the critical slowing 
down for the simulated flat-histogram ensemble and the 
chosen local update dynamics . For the Ising model on 
a square lattice it was shown in Ref. Q that z « 0.74. 

In this paper, we determine the performance of the 
single-spin flip Wang-Landau algorithm for the ferromag- 
netic Ising model as a function of the dimension d. We 
find that the exponent z is a decreasing function of d. 
By using Monte Carlo simulations we determine that 
z « 1.814, 0.74 3, an d 0.438 for d = 1, 2, and 3, respec- 
tively (see Sec. Ill A]) . In the mean- field limit of infinite 
dimension we find that z = using Monte Carlo simu- 
lations and analytic methods (Sec. 1111^ . To round out 
our examination of the single-spin flip algorithm, we also 
estimate z for the q = 10 and q — 20 Potts model in 
d = 2 (see SeclHTl. 

We then discuss two complementary approaches that 
overcome the critical slowing down of the single-spin flip 
flat-histogram algorithm. In Sec. IIIII we demonstrate 
that flat-histogram MC simulations can be speeded up 
by changing the spin dynamics, that is, by using cluster 
updates instead of local updates @, Q • Alternatively, we 
apply the recently developed ensemble optimization tech- 
nique 5] to change the simulated statistical ensemble and 
show that by sampling an optimized histogram instead of 
a flat-histogram, the critical slowing down of the single- 
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spin flip simulation can also be eliminated fSec. lIV|) . For 
the Ising model in one and two dimensions we find that 
both approaches result in optimal A'^^-scaling of the ran- 
dom walk up to logarithmic corrections. 



II. PERFORMANCE OF FLAT-HISTOGRAM 
ALGORITHMS WITH LOCAL UPDATES 

In this section we study the performance of the single- 
spin flip flat-histogram algorithm by measuring the mean 
first-passage time for the d = 1, 2, 3 and mean- field Ising 
models and the q-state Potts model in d = 2. 



A. The mean first-passage times for the Ising 
model in d = 1, 2, and 3 

For the d = 1 Ising model, the mean first-passage time 
time was calculated using the exact density of states and 
50,000 MC measurements for each value of N. A power- 
law fit in the range 70 < iV < 500 gives z = 1.814 ±0.014 
as shown in Fig. ^ However, there is some curvature 
in the data which indicates a deviation from this power- 
law scaling and the effective value of z seems to increase 
as larger system sizes are taken into account. Thus we 
cannot eliminate the possibility that z = 2. 

A crude argument that z < 2 comes from considering 
the two lowest states of a spin chain. The ground state 
has all spins aligned, and the first excited state has a 
single misaligned domain bounded by two domain walls. 
The distance between the domain walls is typically the 
order of L, the linear dimension of the system. For single 
spin flip dynamics the domain walls perform a random 
walk so that the typical number of sweeps required to 
make the transition from the first excited state to the 
ground state scales as L^. If this diffusion time were to 
hold for all energies, the result would he z = 2. However, 
for higher energy states the domain wall diffusion time 
becomes smaller, so z = 2 is an upper bound for the 
d = 1 Ising model. 

The exact density states also was used for d = 2 
The combined results from Ref. Q and this work are 
shown in Fig. |21 where each data point represents a total 
of 50,000 measurements. A power-law fit of the critical 
exponent z for 10 < L < 64 gives z = 0.743 ± 0.002. 

For the d = 3 Ising model the exact Ising density of 
states cannot be calculated exactly, except for very small 
systems. We first used the Wang-Landau algorithm 0,01 
to estimate the density of states. The criteria for the 
completion of each iteration in the determination of g{E) 
was that each energy be visited at least five times and 
that the variance of the histogram be less than 10% of 
the mean. Because g{E) is symmetric about E = 0, 
the calculated values for g{E) and g{—E) were averaged. 
For each value of N ^ g{E) was independently calculated 
ten times and then each g{E) was used to obtain r by 
averaging over 5,000 measurements. The results for the 




FIG. 1: Scaling of the mean first-passage time in the energy 
interval [—N, -\-N] for the flat-histogram random walker of the 
d = 1 Ising model. The dashed line corresponds to a power- 
law fit in the range 70 < A'' < 500, which yields a critical 
exponent oi z = 1.814 ± 0.014. 




FIG. 2: Scaling of the mean first-passage time in the energy 
interval [— 2A'', -|-2A^] for the flat-histogram random walker of 
the d = 2 Ising model. The dashed line corresponds to a 
power-law fit, which yields z = 0.743 ± 0.002. 



first-passage time t were then averaged over ten inde- 
pendent runs. The results are shown in Fig. Inland yield 
z = 0.438 ±0.010. 



B. Mean-Field Ising model 

The results of Sec. Ill Al indicate that the value of z 
decreases with increasing dimension. This dependence 
suggests that it would be interesting to compute r for 
the Ising model in the mean- field limit. We consider the 
"infinite-range" Ising model for which every spin inter- 
acts with every other spin with an interaction strength 
proportional to 1/A^. In this system the energy E and 
magnetization M are simply related, and it is convenient 
to express the density of states in terms of the latter. For 
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FIG. 3: Scaling of the mean first-passage time in the energy 
interval [— SA'^, -|-3Af] for the flat-histogram random walker of 
the d — 3 Ising model. The dashed line corresponds to a 
power-law fit for L > 8, which yields z = 0.743 ± 0.002. 



each value of N we did 50,000 MC measurements of r in 
the range 4 < < 8000. 

We also calculated r from a master equation. Because 
flipping a spin changes the magnetization by ±2, the mas- 
ter equation takes the form (suitably modified near the 
extremes of M = ztN) 

P{M, t + 1) =P(M, t) + T{M + 2, M) P{M + 2, t) 

- [T{M, M + 2)+ T{M, M - 2)] P{M, t) 
+ T{M -2,M)P{M -2,t), (2) 

where P{M, t) is the probability that the system has 
magnetization M at time t and T{Mi, M2) is the prob- 
ability of a transition from a state with magnetization 
Ml to one with M2. The transition probabilities T are 
products of the probability of choosing a spin in the de- 
sired direction times the probability of accepting the flip. 
The latter is the Wang-Landau probability, VF(Mi, Af2), 
which is given by 



(Ml, A/a) = min 



1, 



5(Mi) 



g{M2 



(3) 



The probability of choosing a spin in the desired direction 
is determined as follows. In a transition in which the 
magnetization increases, a down spin must be chosen. 
The probability of choosing a down spin is the number 
of down spins, N^, divided by the total number of spins, 



iV|(M) 
N 



N-M 
2N ■ 



(4a) 



Similarly, in a transition in which the magnetization de- 
creases, an up spin must be chosen. The probability is 



iV|(M) _N + M 



N 



2N 



(4b) 



The transition probabilities can now be written as: 
N^{M + 2) 



T(M + 2, M) 
r(M, M + 2) 
T(M - 2, M) 
r(M, M - 2) 



N 

NjjM) 
N 

Ni{M~2) 
N 

7Vt(M) 
N 



W{M + 2, M) (5a) 

W{M, M + 2) (5b) 

WiM - 2, M) (5c) 

W^(M,M-2) (5d) 



The Wang-Landau probability M^(Mi, M2) can be sim- 
plified for this system in the following way: the density 
of states in terms of the magnetization is 



(6) 



If |M| increases in the transition, that is, |Mi| < IM2I, 
then Vl^(Mi,M2) = 1 because the ratio of the density 
of states exceeds unity. If |M| decreases, there are two 
cases. For M < 0, A^i increases by 1 and M increases by 
2, and we have 



W{M, M -f 2) 



9{M) 



(iv-r+i) 



N + M + 2 



N-M 



(7a) 



g{M + 2) 

If M > 0, decreases by 1 and M decreases by 2, and 
g{M) [nJ N-M + 2 



W{M, M - 2) = 



g{M-2) N + M 



(7b) 

Equations (PJ and Q can be combined into the following 
simple expression 



W{M, M ± 2) = min 



1, 



iV ± M -f 2 



N+\M\ 



(8) 



To compute the mean first-passage time, we take the 
state with magnetization N (the state with the highest 
magnetization) to be absorbing. The initial condition is 



P(Af, 0) 



1 M = -N 
otherwise . 



(9) 



To compute r, we iterate Eq. ^ using Eqs. Q and ^ 
and compute AP{N, t), the change in the value of P{N, t) 
after each iteration t. The mean first-passage time is then 
given by 



T = ^tAP(7V,i). 



(10) 



A comparison of the distribution of first-passage times 
calculated from the master equation and from a direct 
simulation is shown in Fig. 2| As shown in Fig. |S1 the 
scaling of the mean first-passage time with N is consis- 
tent with t/N'^ ~ In TV. The fit includes both the master 
equation and simulation results. 
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FIG. 4: Comparison of the distribution of first-passage times 
for the mean-field Ising model computed by iterating the 
master equation (line) and by direct simulation (points) for 
— 32. The master equation results were normalized to 
50,000 MC estimates. 



The logarithmic behavior of t can be determined ana- 
lytically from the transition rates using the first-passage 
time methods described in Ref. Let i(m) be the 

mean first-passage time starting from magnetization m = 
M/N and ending at the absorbing boundary m = 1 with 
a reflecting boundary at m = —1. We seek r !)• 
The mean first-passage time satisfies, 

tlm) = [t{m — Sm) + 5t\T^{m) 

+ [tlm + 5m) + 5t]T^{m) 

+ [t{m)+5t]To{m), (11) 

where r_(m) = T{M,M - 2), r+(m) = T{M,M + 2) 
are the transition probabilities, Tq = 1 — r_ — T+ is the 
waiting probability, Sm ~ 2/N is the magnetization step 
size, and = 1 is the time unit for a step. Equation 1)11(1 
can be cast in differential form by expanding t{m) to 
second order in Sm, 



where 



1 + h{m)t'{m) + f{m)t"{m) = 0, (12) 



/(m) = (W/25i)[T+(TO) +T_(m)], (13) 



and h{m) — {Sm/ 5t)[T^{m) — T_(m)]. The absorbing 
boundary requires that 



i(l) = 0. 



(14) 



At the reflecting boundary, the transition to the left van- 
ishes, T_(— 1) = 0, and we have 

t{-l) ^ [t{-l + Sm) + St]T+{-l) 

+ [t{-l) + St][l^T+i^l)]. (15) 

If we expand Eq. (|15|) to first order in Sm, we obtain the 
boundary condition 



t'{-l) 



2/(-l)' 



(16) 



where /(-I) = {Sm"^ /2St)T+{m) 

The flat-histogram feature of the Wang-Landau 
method means that the random walk satisfies the sim- 
ple detailed balance condition. 



T+ (m) — T_ (m + Sm) 
T^{m) = T^{m^Sm). 



(17a) 
(17b) 



We expand the detailed balance conditions to first or- 
der in Sm and find a relation between the coefficients in 
Eq. 113, 

h{m) = f{m), (18) 

with the result that Eq. (|12|l reduces to 

1 + f'{m)t'{m) + f{m)t"{m) =0. (19) 

If we take into account the reflecting boundary condition, 
the solution to the differential equation is 



x + l + Srn/2 

t(m) = T - — dx, 

1-1 fix) 



(20) 



where t = t{—l) is the first-passage time we are seeking. 
We use the absorbing boundary condition to find that 



^ X + 1 + Sm/2 



dx. 



(21) 



Because /(m) is an even function, the first term in the 
integrand vanishes. We substitute Sm — 2/N and use 
Eq. H13|) for /(m) to obtain the leading behavior of r. 



r = 2N^ 



-1 1 - n{x) 



-dx. 



(22) 



From Eqs. Q, I^Jl, and (jSJl we have that To{m) = \m\ 
1/N, which yields the desired result 



r = 2N' 



dm 



l-\m\ + 1/N 



N^logN. (23) 



Note that the integral can be interpreted as the average 
waiting time and that this average is dominated by the 
two boundaries. 



C. The Potts model 

In Ref. it was shown that the non-vanishing scal- 
ing exponent z quantifies the critical slowing down of 
the flat-histogram simulation near a critical point. Most 
strikingly, this slowing down can be observed as a sup- 
pression of the local diffusivity of the flat-histogram ran- 
dom walker (see Sec. IIV|) . The Ising model undergoes a 
second order phase transition. Here we extend our analy- 
sis to flrst-order phase transitions as is observed in g-state 
Potts models with g > 4 on a square lattice. 
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FIG. 5: Scaling of the mean first-passage time in the magneti- 
zation interval [— Af, -\-N\ for the flat-histogram random walk 
of the mean-field Ising model. Note the linear scale on the 
y-axis. The -I- symbol denotes the simulation results and x 
denotes results derived from the master equation. The dashed 
line corresponds to a logarithmic fit of the combined results 
for L > 50. The slope is found to be 0.977 ± 0.004. 



FIG. 6: Scaling of the mean first-passage time in the energy 
interval [— 27V, 0] for the flat-histogram random walk of the 
d = 2, g = 10 and g = 20 Potts models. The dashed lines 
correspond to power-law flts for L > 8 and yield z = 0.786 ± 
0.010 and z = 0.961 ± 0.008 for the g = 10 and g = 20 Potts 
models, respectively. 



Because the exact density of states is not known for 
the g-state Potts model, we follow the same procedure as 
for the Ising model in d = 3 making use of the Wang- 
Landau algorithm. Based on these estimates for the den- 
sity of states, we measured the mean first-passage time 
as a function of the hnear system size L as is illustrated 
in Fig. El For q — 10 we find that the critical exponent 
z becomes z ~ 0.786 ± 0.010, which is close to our re- 
sult for the d — 2 Ising model thereby showing the close 
proximity of this weak first-order phase transition to a 
continuous phase transition. For the g = 20 Potts model 
we find z = 0.961 ± 0.008, illustrating increasing critical 
slowing down as the strength of the first-order transition 
increases. However, we cannot rule out exponential slow- 
ing down for even larger values of q as was observed for 
multicanonical simulations of droplet condensation in the 
d = 2 Ising model i2j ■ 



III. CHANGING THE DYNAMICS: 
CLUSTER UPDATES 

In Sec. ^ we found that single-spin flip fiat-histogram 
Monte Carlo simulations suffer from a critical slowing 
down, that is, z > 0, for the Ising model in d = 1, 2, 
and 3. The critical slowing down associated with single- 
spin flip algorithms near critical points in the canonical 
ensemble has been reduced by the introduction of clus- 
ter algorithms such as the Swendsen-Wang algorithm 
and the Wolff algorithm In this section we follow 

a similar approach and combine efficient cluster updates 
with flat-histogram simulations. 

In conventional cluster algorithms |^ ^3 simulating 
the canonical ensemble, clusters of parallel spins are built 



up by adding aligned spins to the cluster with a proba- 
bility p(/3), which explicitly depends on the temperature 
/3 = 1/T. When simulating a broad- histogram ensem- 
ble in energy there is no explicit notion of temperature, 
and thus no straightforward analogue of a cluster algo- 
rithm. Recently, Reynal and Diep suggested a solution of 
this problem by using an estimate of the microcanonical 
temperature f3{E) = dS{E)/dE 

A complementary approach is to sample a broad- 
histogram ensemble using an alternative representation 
of the system's partition function for which it is possible 
to genuinely introduce cluster updates. An example is 
the multibondic method introduced by Janke and Kap- 
pler Q which uses the Fortuin-Kasteleyn (FK) represen- 
tation [m in the context of multicanonical sampling. Al- 
though the multibondic method performs local updates of 
the graph in the FK representation and cluster updates 
of the spin configurations, Yamaguchi and Kawashima 
were able to show that a global and rejection free update 
of the graph in the FK representation is possible 7] . We 
will discuss methods based on graph representations in 

sec inra 

The above methods can only be applied to Ising and 
Potts models. In Sec HTTBl we introduce a new represen- 
tation using multigraphs that allows cluster updates for 
continuous spin models in the context of broad- histogram 
ensembles. 

In the following we assume that the Hamiltonian of the 
system has the form 

H[<J] = -Y,h^JX^l]^ (24) 

where cr denotes the configuration of the system, Ui the 
spin of site i, and the sum is over all bonds {ij) that 
define the lattice of the system. 
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A. Graph representation 

We first review the cluster methods based on an addi- 
tional graph variable and present results showing that 
z = for multibondic flat-histogram simulations of 
Ising models. In the spin representation the canonical 
partition function of the Ising-Potts model is given by 
Z = X^CT l^crCc), where the summation is over all spin 
configurations, and the weight of a spin configuration is 

The g-state Ising-Potts system also can be represented 
by a sum over all graphs lo that can be embedded into 
the lattice. The partition function is then 

Z = Y,WU^), (25) 

and the weight of a graph uj is given by 

T4^^(u;) (26) 

n{uj) is the number of of bonds in the graph lo, c{ll!) 
the number of connected components (counting isolated 
sites), and rtf, is the total number of bonds in the lattice. 
Because the graph uj can be viewed as a subgraph of the 
lattice, its bonds are sometimes referred to as "occupied 
bonds," and the bond probability of the ferromagnetic 
Ising-Potts model with hij[cri,crj] = Sai,crj is given byna 

p(/3) = 1 - (27) 

A third representation of the Ising-Potts model is 
the spin-bond representation, which is employed in the 
Swendsen-Wang algorithm |9|. In the spin-bond repre- 
sentation the system is characterized by both spins and 
bonds, with the requirement that a bond can be occupied 
only if it is satisfied, that is if the two spins connected by 
a bond in oj have the same value. In this representation 
the partition function is 

^ = 5]W^aa,('T,Cc.), (28) 

and the weight is given by 

W„^{a,uj)^p"^^\l-pr-"^^^A{a,Lj), (29) 

where A(ct, w) = 1 if all occupied bonds are satisfied and 
zero otherwise. 

It is natural to introduce a density of states for each 
of the three representations: For the spin representation 
in terms of the energy E, we write 

gAE)=Y,mW]-E). (30) 

For the bond representation in terms of the number of 
occupied bonds n and number of clusters c, we have 

g^{n, c) = ^ 6{n{uj) - n)5{c{uj) - c). (31) 




6-1 1 1 1 1 r- 

32 64 128 256 512 



FIG. 7: Scaling of the mean first-passage time applying clus- 
ter updates in the bond representation for the d = 1 Ising 
model. Note the linear scale of the y-axis. 

And for the spin-bond representation in terms of the 
number of occupied bonds n, we write 

5ac.(n) = ^5(n(c^)-n)A(a,cc;). (32) 

UJ,(T 

The corresponding forms of the partition function are 
Z = ^5,(i?)e-''^ (spin) (33) 

E 

^ = E5-Kc)p"(l-p)"^-"g= (bond) (34) 

n,c 

^ = ^5<xc.(n)p"(l-Pr''~"- (spin-bond) (35) 

n 

The Wang-Landau algorithm can be applied in all three 
representations. As discussed in Sec. the spin repre- 
sentation has a relatively large value of z. 

The bond representation generally requires a two- 
dimensional histogram, but for the d — 1 Ising-Potts 
model, the number of clusters c is completely deter- 
mined by the number of occupied bonds n and a one- 
dimensional histogram is sufficient. We simulated the 

= 1 Ising model in the bond representation and found 
that the mean first-passage time scales as TV^logL, as 
shown in Fig. |7| It is noteworthy that the domain wall 
arguments used to explain the large value of z in one 
dimension do not apply in the bond representation. 

In higher dimensions n does not determine c, and it 
is simpler to use the spin-bond representation. The re- 
sulting multibondic algorithm Q is very similar to the 
Swendsen-Wang algorithm, and can be summarized as 
follows, 

1. Choose a bond at random. 

2. If the bond is satisfied and occupied (unoccupied), 
make it unoccupied (occupied) with probability 
min[gi^a-in)/guja-in'),l], where n and n' are the 
number of occupied bonds before and after the 
change, respectively. Update the observables. 
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3. If the bond is unsatisfied, the system stays in its 
original state and the observables are updated. 

4. After one sweep of the lattice, identify clusters 
(connected components) and assign spin values to 
them with equal probability. 

The algorithm is ergodic and satisfies detailed bal- 
ance. The ergodicity is obvious from the fact that the 
Swendsen-Wang algorithm is ergodic. In general, the de- 
tailed balance relation can be written as, 



Tr{a,uj)p{a a'\uj)p{uj ^'\<7') 
= Tr{a' ,Lu')p{uj' — > uj\a')p{a' — > cr\Lu), 



(36) 



where 7r(cr, u) is the equilibrium probability of microstate 
(cr, w), and p(cr —>■ cr'\uj) is the transition probability from 
a to a' given that the bond configuration is fixed at to. 
The algorithm is designed to sample the distribution, 



7r(cr,tj) 



1 



The transition probability of a spin flip is 
„ . 1 



A(a,t^). 



p{a 



A(a,a;) 



(37) 



(38) 



and the transition probability of a bond change (lo' ^ lS) 
is given by. 



p{u} 



1 



gau>{n{uj)) 
gau>{n{u}')) 



(39) 



because the probability of choosing a bond to change 
is l/rib, and the probability of accepting the change is 
min [(7ct(.j('t-(w))/(7ct(.j("-('^'))7 l] • It is easy to verify that 
the detailed balance relation, Eq. is satisfied. 

We have applied this algorithm to the d = 2 Ising 
model and measured the mean first-passage time, which 
also scales as iV^logL, as shown in Fig.|S| At least 10^ 
first-passage times were measured for each value of L. 



B. Spin-multigraph representation 

We now consider a representation of the partition 
function that can be used to implement cluster updates 
in fiat-histogram simulations of continuous spin models 
such as 0(n) models. Like the spin-bond representation 
it has a configuration and a graph variable. In contrast to 
the spin-bond representation, the graph is a multigraph 
and can therefore have multiple bonds between two ver- 
tices. Before introducing this representation, we briefly 
review how a simple high-temperature series representa- 
tion can be sampled in a flat-histogram simulation. 

We write the partition function as a series in the inverse 
temperature P 



(40) 




FIG. 8: Scaling of the mean first-passage time using cluster 
updates in the spin-bond representation of the d = 2 Ising 
model. Note the linear scale of the y-axis. 



with a corresponding density of states given by 



9P\ 



(41) 



If H[a\ < 0, then the weight of a conflguration-order pair 



(42) 



is always positive, and we can perform a flat-histogram 
simulation in the extended phase space (cr, n) using a 
and n updates. If the a and n updates are performed 
independently of each other, an n update n ^ n' (for 
example n' = n ± I) is accepted with probability 



Pr, 



n'\gf}[n') 



(43) 



A configuration update ct ^ cr' (for example the change 
of the configuration of a single site Ci — > cr^) at fixed n is 
accepted with probability 



P. 



H[a] 



1 



(44) 



In a practical computation n is truncated at some cutoff 
A restricting the largest inverse temperature that can be 
reached to /3max ^ A/T^ T^, where V is the volume of 
the system. 

Because the weight of a configuration as defined in 
Eq. (|42|l is not a product of bond weights on the lat- 
tice, the (cr, n)-representation cannot be directly used to 
implement cluster updates. To obtain such a representa- 
tion we start from 



(45) 



8 



and expand each exponential term in /3 with a separate 
integer variable tOij 



(46) 



We identify each set of ojij on the lattice with a multi- 
graph cu that has iVij bonds between site i and site j. If 
we write n{u)) = ^ Uij for the total number of bonds in 
the multigraph w, gp{n) as defined in Eq. H41|l is given 

by 
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where the sum is over all spin configurations a and all 
multigraphs uj on the lattice. The corresponding parti- 
tion function is given by Eq. (|40|l . and the weight of a 
configuration-multigraph pair is given by 



0.5 



FIG. 9: Scaling of the mean first-passage time for the d — 2 
Ising model using the simple high temperature representation 
with single spin updates (•) and using the spin-multigraph 
representation with Swendsen-Wang cluster updates (□). The 
lines are drawn as a guide to the eye. 



IT";3 (O-, W) = Wij {(Ji , (Tj , LUij 



(48) 



either none or both sites are fiipped and Wij the weight 
after a flip of one of the sites (note that this requires 
Wtj {<J^,<7j) = W,j {a[ , ct' ) and W^j {a[,<jj) = Wij (tXi , tr' ) ) , 



where Wij{ai,aj,ujij) = hij[ai,aj]'^^^ /ujijl is the weight 

of the bond (ij) in the lattice. The energies in the sys- ^^^^ probability of adding a site to^a'cluster is giVen by 
tem have to be shifted so that hij\ai,ai] > to ensure 



that all weights are positive. Again two types of updates, 
one in the multigraph uj and one in the conflguration a 
are needed. Before discussing them in detail, we men- 
tion some differences of the spin-multigraph representa- 
tion compared to the representations used in Sec. IIII Al 
In the spin-multigraph representation we cannot in gen- 
eral integrate explicitly over all configurations to obtain 
a multigraph-only representation similar to the bond rep- 
resentation except for Ising or Potts models. It also is not 
possible to apply the global graph update of Yamaguchi 
and Kawashima to the multigraph so that local up- 
dates have to be used. Finally the representation can be 
viewed as the classical limit of the stochastic series expan- 
sion method (SSE) for quantum systems [l7( . The order 
of the operators in the operator-string used in the SSE is 
unimportant for a classical system, so that one only has 
to count how often an operator hij[ai,aj] occurs (w.y). 

If we use a fiat-histogram algorithm in the order n, a 
multigraph update uj ^ uj' (as before uj'^j — ujij ± 1 for 
one randomly chosen bond (ij) is a good choice) will be 
accepted with probability 



= 1 — min 



-,1 



(50) 



For the canonical ensemble simulation of the Ising model 
with Ty(t,T,/9) = and VF(T,i,/3) = 1, Eq. ^ re- 
duces to the well known \ — for parallel spins as 
discussed before. In the spin-multigraph representation 
the probability to add a site to a cluster is 



1 



(51) 



Thus only sites that are connected by at least one bond 
of the multigraph can be added to the cluster. For the 
Ising model with hij[(7i, Uj] — Sa-^.a^ , the weight of a con- 
figuration (cr, uj) as defined in Eq. (|48|l is 



Wii{a,uj) = A(cr,tj) 



n 



(52) 



gpjnjuj)) yr Wtjiat,aj,uj'ij) ^ 



(49) 



As the weight of a spin-multigraph configuration is a 
product of bond weights, a cluster update scheme such 
as the Swendsen-Wang Q or Wolff algorithm can be 
used to update the spin configuration. After choosing 
a flip operation for the spins, the only difference to the 
canonical cluster algorithms is the probability of adding 
a site to a Swendsen-Wang cluster: If Wij is the weight if 



where A(cr, cj) — 1 if all bonds in the multigraph uj are 
satisfied and otherwise, and Pconnoct simplifies to 1 if 
two sites are connected by one or more bonds in the 
multigraph and to if they are not. 0{n) models can 
be simulated by performing the spin update with respect 
to a mirror plane that is randomly chosen for one cluster 
update |T(tI |. 

In a local update scheme, a configuration update cr — > 
cr' at fixed uj will be accepted with probability 



{ij) 



hij[ai,aj] 



,1 



(53) 
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FIG. 10: Scaling of the mean first-passage time to reach a 
low energy state for the d = 1 XY model using single spin 
updates (•) and Swendsen-Wang cluster updates (□) in the 
spin-multigraph. Lines are drawn as a guide to the eye. 




FIG. 11; Scaling of the mean first-passage time of the opti- 
mized ensemble for the d=\ Ising model using Metropolis up- 
dates. The solid line corresponds to a power-law fit with expo- 
nent 0.972 ± 0.002. The inset shows the optimized histogram 
which exhibits a power-law divergence — _Eo)/A'']^'^°. 



We have applied a flat-histogram sampling of this rep- 
resentation to the d — 2 Ising model with periodic bound- 
ary conditions and to the d ~ 1 XY model with open 
boundary conditions. To compare the effectiveness of 
the cluster updates to the single-spin flip updates, we 
measure the mean first-passage time Tspin from a random 
state cr at n = 0, which is equivalent to = for all 
(ij), to an ordered state. For the Ising model this state 
is chosen to be the ground state and for the d = 1 model 
any state which is within 2.5% of the width of the spec- 
trum to the ground state. We perform a fixed ratio of 
spin to multigraph updates and measure the mean first- 
passage time in spin updates for different ratios of spin 
to multigraph updates, rgpin converges to a fixed number 
with an increasing fraction of multigraph updates, and 
for sufficiently many multigraph updates we can assume 
the multigraph is equilibrated for a particular configura- 
tion. The cutoff A is chosen so that the average n for 
which the ground state is reached is much smaller than 
A. 

Figure shows the scaling of Tspm for the d = 2 Ising 
model for single-spin flips in the simple representation 
of Eq. (|42|l and for Swendsen-Wang cluster updates in 
the spin-multigraph representation. The single-spin flip 
updates show a power law behavior with z = 0.85 ± 0.06, 
while the cluster updates can be fitted by a logarithmic 
behavior or a power law with z = 0.10 ± 0.02. Figure [TUI 
shows the scaling of Tgpin for the d = \ XY model using 
single spin updates and Swendsen-Wang cluster updates 
in the spin-multigraph representation. 



IV. CHANGING THE ENSEMBLE: 
OPTIMIZING THE SAMPLED HISTOGRAM 

An alternative to changing the dynamics from local to 
cluster updates for flat- histogram sampling is to optimize 



the simulated statistical ensemble and retain local spin- 
flip updates. When simulating a flat-histogram ensem- 
ble, the random walker is slowed down close to a critical 
point, which is reflected in the suppressed local diffusiv- 
ity. To overcome this critical slowing down, we can use 
this information and define a new statistical ensemble 
w{E) by feeding back the local diffusivity D{E) 5j. Af- 
ter the feedback additional weight is shifted toward the 
critical point. The new histogram is no longer flat, but 
exhibits a peak at the critical energy, that is, resources 
in the form of local updates are shifted toward the criti- 
cal point. Ultimately, this feedback procedure results in 
an optimal histogram n„(i5) that is proportional to the 
inverse of the square root of the local diffusivity D{E) 



1 



VD{E) 



(54) 



For the d = 2 Ising model it was shown in Ref. 
that the mean first-passage time scales as r oc (iVlnA^)^ 
for the optimized ensemble and satisfies the scaling of 
an unbiased Markovian random walk up to logarith- 
mic corrections. It also was found that the distribution 
of statistical errors in the calculated density of states 
g{E) = nw{E)/w{E) is uniformly distributed in energy, 
in contrast to calculations based on flat-histogram sam- 
pling. 

Here we show results for the optimized ensemble of the 
d = 1 Ising model for both Metropolis and A^-fold way 
updates For single-spin flip Metropolis updates, we 
find that the histogram is shifted toward the ground-state 
energy Eq = —N and follows a power-law divergence, 
n^{E) cx [{E-Eo)/N]^-^°'^° °^ (see the inset of Fig.HIJ- 
However, the mean first-passage time of the random walk 
in the energy interval [—N, 0] still exhibits a power-law 
slowdown with z = 0.972 ± 0.002 as illustrated in the 
main panel of Fig. 1111 This remaining slowdown may 
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FIG. 12: Scaling of the mean first-passage time of the op- 
timized ensemble for the ID Ising model using Ai'-fold way 
updates (note the linear scale on the y-axis). The solid line 
corresponds to a logarithmic fit. The inset shows the opti- 
mized histogram. 

originate from the slow dynamics that occurs whenever 
two domain walls reside on neighboring bonds, and a spin 
flip of the intermediate spin is suggested with a proba- 
bility of 1/A^. This argument suggests that changing the 
dynamics from simple Metropolis updates to TV-fold way 
updates would strongly increase the probability of anni- 
hilating two neighboring domain walls. 

We optimized the ensemble for TV-fold way updates and 
found that the mean first-passage time is further reduced 
and scales as t c>c (iVlnA'^)^ as shown in Fig. ^1 This 
scaling behavior corresponds to the results for the opti- 
mized ensemble for the d = 2 Ising model Q. However, 
for the d ^ 2 model the optimized ensembles for both 
Metropolis and iV-fold way updates resulted in the same 
scaling behavior of the form t oc (A^ IniV)^. 



V. CONCLUSIONS 

We have studied the performance of flat-histogram 
simulations of Ising and Potts models in different dimen- 
sions. For one, two, and three dimensions the simulations 
show critical slowing down with a critical exponent z > 0. 
This behavior is analogous to the critical slowing down 
of canonical ensemble simulations of these models using 
single spin updates. Our numerical results show that z 
decreases as a function of the dimensionality d and van- 
ishes (up to logarithmic corrections) for the infinite range 
Ising model. 

We demonstrated that the critical slowing down of 
the flat-histogram simulations can be overcome by ei- 
ther changing the representation of the system to allow 
cluster updates or by changing the simulated ensemble 
and keeping local updates. In order to apply cluster 
updates to continuous spin models in connection with 
broad-histogram simulations, we introduced a new spin- 
multigraph based representation. The broad-histogram 
simulations which take advantage of cluster updates or 
an optimized statistical ensemble do not suffer from criti- 
cal slowing down and show optimal scaling in comparison 
to an unbiased Markovian random walk. 
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